A personalized and non-parametric framework for detecting changes in gait cycles

Department of Mechanical & Aerospace Engineering, University at Buffalo

Department of Mechanical & Aerospace Engineering, University at Buffalo

Department of Industrial and Systems Engineering, University at Buffalo

Farmer School of Business, Miami University

Farmer School of Business, Miami University

---

1 R Setup and Required Packages

In the following code chunk we load the packages used to support our analysis.

# check if packages are not installed; if yes, install missing packages
packages = c("tidyverse", "magrittr", # typical data analysis packages
             "MALDIquant", # match closest points between two vectors
             "foreach", "parallel", # packages used for parallelizing some code chunks
             "R.matlab")
newPackages = packages[!(packages %in% installed.packages()[,"Package"])]
if(length(newPackages) > 0) install.packages(newPackages)

# using the library command to load all packages; invisible used to avoid printing all packages and dependencies used
invisible(lapply(packages, library, character.only = TRUE))

# source("./Functions.R") # our custom built functions

# set.seed(2020)
# startTime <- Sys.time()

---

2 Gait Acceleration Data and Rational Subgroups

We converted the raw IMU acceleration signals in the local to the global reference frame and removed the gravity. The acceleration signals were then transformed back to the local reference frame. Sagittal acceleration, lateral acceleration, and acceleration magnitude signals were then calculated. We used the vertical acceleration component to segment the gait cycles in order to isolate individual gait cycles. The acceleration profiles along with the experimental time stamps were stored in mat files.

In this code chunk we load the segmented acceleration mat files. The \(1^{st}\) 10 minutes of the data were considered as warm up period and thus excluded form the analysis. The acceleration profiles of the gait cycles during each walking cycle were grouped into a rational subgroup. We stored the start and end times of the walking cycles in csv files and loaded and used them in the following code chunk. The acceleration profiles within each start and end times to subgroup the gait cycles.

for (id in setdiff(1:15, 13)) {
  #################### Read start and end time of the subgroups
  video <- read.csv(paste0(file="../Data/csvFiles/Sub", id, ".csv"))

  video_leave <- video$Leaves
  video_leave <- video_leave[!is.na(video_leave)]
  video_leave <- video_leave - 600
  video_leave <- video_leave[video_leave > 0]

  video_enter <- video$Enters
  video_enter <- video_enter[!is.na(video_enter)]
  video_enter <- video_enter - 600
  video_enter <- video_enter[video_enter > 0]

  ################################################### Read from mat files
  raw1 <- readMat(paste0("../Data/matFiles/Subject", id, "_aZ_seg.mat"))

  raw2 <- raw1$gait

  num_rows <- length(raw2)/5

  aM <- list("vector")
  aS <- list("vector")
  aL <- list("vector")
  exp_time <- c()
  for (i in 1:num_rows) {
    aM[[i]] <- raw2[[1*num_rows + i]][[1]][1,]
    aS[[i]] <- raw2[[2*num_rows + i]][[1]][1,]
    aL[[i]] <- raw2[[3*num_rows + i]][[1]][1,]
    exp_time[i] <- raw2[[i+4*num_rows]][[1]][1,1]
  }

  ###################################################### Create subgroups
  if (video_leave[1] > video_enter[1]){
    video_leave <- c(exp_time[1], video_leave)
  }

  if (tail(video_leave, 1) > tail(video_enter, 1)){
    video_enter <- c(video_enter, tail(exp_time, 1))
  }

  conf <- 0.05 * mean(video_enter - video_leave) # 5% confidence for each subgroup

  # Match video subroup times against the experimental times from mat files
  leave_index <- match.closest(video_leave + conf, exp_time)
  enter_index <- match.closest(video_enter - conf, exp_time)

  aM_list <- list()
  aS_list <- list()
  aL_list <- list()
  t_list <- list()
  for (i in 1:length(leave_index)) {
    aM_list[[i]] <- aM[leave_index[i]:enter_index[i]]
    aS_list[[i]] <- aS[leave_index[i]:enter_index[i]]
    aL_list[[i]] <- aL[leave_index[i]:enter_index[i]]
    t_list[[i]] <- exp_time[leave_index[i]:enter_index[i]]
  }

  len <- c()
  for (i in 1:length(aM_list)) {
    len <- c(len, lengths(aM_list[[i]]))
  }

  cut_len <- quantile(len, 0.02)

  for (i in 1:length(aM_list)) {
    remove_ind <- which(lengths(aM_list[[i]]) < cut_len)
    remove_ind <- c(remove_ind, 3000)
    aM_list[[i]] <- aM_list[[i]][-remove_ind]
    aS_list[[i]] <- aS_list[[i]][-remove_ind]
    aL_list[[i]] <- aL_list[[i]][-remove_ind]
    t_list[[i]] <- t_list[[i]][-remove_ind]
  }

  batchSize <- lengths(aM_list)

  ###############################################
  aM_list_all <- do.call(c, aM_list)
  aS_list_all <- do.call(c, aS_list)
  aL_list_all <- do.call(c, aL_list)
  t_all <- do.call(c, t_list)
  ###############################################
  aM_mat <- do.call(rbind, aM_list_all)[,1:cut_len]
  aS_mat <- do.call(rbind, aS_list_all)[,1:cut_len]
  aL_mat <- do.call(rbind, aL_list_all)[,1:cut_len]

  ##################
  assign(paste0("sub", id, "_aMSLT"),
         list(aMag=aM_mat, aSag=aS_mat, aLat=aL_mat, time=t_all, batchSize=batchSize))
}

########################################### Save
save(sub1_aMSLT, sub2_aMSLT, sub3_aMSLT, sub4_aMSLT, sub5_aMSLT, sub6_aMSLT,
     sub7_aMSLT, sub8_aMSLT, sub9_aMSLT, sub10_aMSLT, sub11_aMSLT, sub12_aMSLT,
     sub14_aMSLT, sub15_aMSLT,

     file="../Data/rData/subGs_for_depth.Rdata")

3 Depth Calculation

Explain about depth concept in general; centrality and outlyingness

3.1 Mode Depth

Calculate mode depth…

load(file="../Data/rData/subGs_for_depth.Rdata")

cores = detectCores() - 2 # to give the server some breathing Room
cl = makePSOCKcluster(cores)
registerDoParallel(cl)

for (id in setdiff(1:15, 13)) {
  
  aMag <- get(paste0("sub", id, "_aMSLT"))$aMag
  tExp <- get(paste0("sub", id, "_aMSLT"))$time
  
  ################################## Mode Depth
  inControl <- aMag[1:500,]
  inControlDepth <- {depth.mode(fdata(inControl))}$dep
  
  ###################### foreach
  onlineDepth <- c()
  end <- nrow(aMag)
  onlineDepth <- foreach(i=501:end, .packages = c('fda.usc', 'tidyverse'),
                         .combine='c') %dopar% {
                           temp1 <- aMag %>% .[i,]
                           append1 <- rbind(inControl, temp1)
                           temp2 <- depth.mode(fdata(append1))
                           
                           temp3 <- temp2$dep %>% .[501]
                           temp3
                         }
  
  modeDepth <- c(inControlDepth, onlineDepth)
  
  assign(paste0("sub", id, "_mode_mag"),
         list(modeDepth=modeDepth, tExp=tExp))
}

save(sub1_mode_mag, sub2_mode_mag, sub3_mode_mag, sub4_mode_mag, sub5_mode_mag, sub6_mode_mag,
     sub7_mode_mag, sub8_mode_mag, sub9_mode_mag, sub10_mode_mag, sub11_mode_mag, sub12_mode_mag,
     sub14_mode_mag, sub15_mode_mag,
     file="../Data/rData/mode_mag.Rdata")

3.2 MFHD Depth

Calculate MFHD depth…

Explain since it take very long it was performed on supercomputer cluster in parallel and the eval option in the following code chunk was set to false since this could take weeks to run on the local computer.

#################################################
load(file="../Data/rData/subGs_for_depth.Rdata")

cores = detectCores() - 2 # to give the server some breathing Room
cl = makePSOCKcluster(cores)
registerDoParallel(cl)

for (id in setdiff(1:15, 13)) {
  
  aSag <- get(paste0("sub", id, "_aMSLT"))$aSag
  aLat <- get(paste0("sub", id, "_aMSLT"))$aLat
  tExp <- get(paste0("sub", id, "_aMSLT"))$time
  
  ################################## Mode Depth
  inControlSag <- aSag[1:500,]
  inControlLat <- aLat[1:500,]
  inControlDepth <- {MFHD(y1=inControlSag, y2=inControlLat, alpha=0.125, Beta=0.5)}$MFHDdepth[1,]
  
  ###################### foreach
  onlineDepth <- c()
  end <- nrow(aSag)
  onlineDepth <- foreach(i=501:end, .packages = c('fda.usc', 'tidyverse', 'MFHD'),
                         .combine='c') %dopar% {
                           temp1 <- aSag %>% .[i,]
                           temp2 <- aLat %>% .[i,]
                           
                           append1 <- rbind(inControlSag, temp1)
                           append2 <- rbind(inControlLat, temp2)
                           
                           temp3 <- MFHD(y1=append1, y2=append2, alpha=0.125, Beta=0.5)
                           
                           temp4 <- temp3$MFHDdepth %>% .[1,501]
                           temp4
                         }
  
  MFHDdepth <- c(inControlDepth, onlineDepth)
  
  assign(paste0("sub", id, "_MFHD"),
         list(MFHDdepth=MFHDdepth, tExp=tExp))
}

save(sub1_MFHD, sub2_MFHD, sub3_MFHD, sub4_MFHD, sub5_MFHD, sub6_MFHD, sub7_MFHD, sub8_MFHD,
     sub9_MFHD, sub10_MFHD, sub11_MFHD, sub12_MFHD, sub14_MFHD, sub15_MFHD,
     file="../Data/rData/MFHD.Rdata")

3.3 Visualize the MFHD Depth

load(file="../Data/rData/MFHD.Rdata")

for (id in setdiff(1:15, 13)) {
  dep <- get(paste0("sub", id, "_MFHD"))$MFHDdepth
  t <- get(paste0("sub", id, "_MFHD"))$tExp
  
  cat("###", paste0("Subject", id), "{-}",'\n')
  plot(t, dep, pch=16, cex=1,col=c(rep("red", 500), rep("black", length(dep)-500)))
  legend("topright", legend=c("Baseline data", "New data"), col=c("red", "black"),
         pch=16)
  cat('\n \n')
  
}

Subject1

Subject2

Subject3

Subject4

Subject5

Subject6

Subject7

Subject8

Subject9

Subject10

Subject11

Subject12

Subject14

Subject15

---

References

---

1. Email: saebraga@buffalo.edu | Website: LinkedIn↩︎

2. Email: jiyeonk@buffalo.edu | Phone: +1-716-645-6063 | Website: University at Buffalo Official↩︎

3. Email: loracavu@buffalo.edu | Phone: +1-716-645-4696 | Website: University at Buffalo Official↩︎

4. Email: fmegahed@miamioh.edu | Phone: +1-513-529-4185 | Website: Miami University Official↩︎

5. Email: farmerl2@miamioh.edu | Phone: +1-513-529-4823 | Website: Miami University Official↩︎